if (Cobb_Douglas)
    tech='_Cobb_Douglas';    
    load ../temp/share_model_age.txt
    share_model_age=reshape(share_model_age, [2 2 T 3]);
else
    tech='';            
end

if (quo)    
    pbar_age=zeros(N,T);    
else
    
    insample_baseline=load(['../temp/insample_quo' tech '.txt']);
    pbar_age=load(['../temp/pbar_age' tech '.txt']);
    
    % prices before the price change
    pbar_age=reshape(pbar_age, [N T]);
    
    if (reduce_W)
        name='W';
        W=(1-pct_reduction/100)*W;
        W_m=(1-pct_reduction/100)*W_m;
        W_f=(1-pct_reduction/100)*W_f;
    elseif (reduce_W_constant_full_income)
        name='W_constant_full_income';        
        W_m=(1-pct_reduction/100)*W_m;
        W_f=(1-pct_reduction/100)*W_f;
    elseif (reduce_P_c)
        name='P_c';        
        P_c=(1-pct_reduction/100)*P_c;
    elseif (reduce_p)
        name='p';        
        p=(1-pct_reduction/100)*p;
    end
end

% reset share parameters
age(1:N)=5;
nkid(1:N)=1;
nkid_0_5(1:N)=1;

set_regressors;
evaluate;

if (quo)
    dlmwrite(['../temp/insample_quo' tech '.txt'],insample);    
    insample_baseline=insample;
end

%% outcomes
outcome=zeros(11,2);

for m=1:2
    j=0;
    mp=m-1;

    % expenditure
    j=j+1;
    eval(['outcome(j,m)=mean(pbar_X(married==' int2str(mp) '),''omitnan'');']);

    % price
    j=j+1;
    eval(['outcome(j,m)=mean(pbar(married==' int2str(mp) '),''omitnan'');']);
    
    % quantity (5)
    j=j+1;    
    eval(['outcome(j,m)=mean(tau_m(married==' int2str(mp) '),''omitnan'');']);

    j=j+1;
    eval(['outcome(j,m)=mean(tau_f(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1;
    eval(['outcome(j,m)=mean(g(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1;
    eval(['outcome(j,m)=mean(Y_c(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1;
    eval(['outcome(j,m)=mean(X(married==' int2str(mp) '),''omitnan'');']);    

    % so far there are 7 outcomes
        
    % log achievement at age 6 (age 5 investment)
    j=j+1; % 8
    eval(['outcome(j,m)=mean(qd_logPsi(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1; % consumption unit (9)
    aux=alpha.*qd_logPsi;
    eval(['outcome(j,m)=mean(aux(married==' int2str(mp) '),''omitnan'');']);
    
    max_j=j;
end    

logPsi=zeros(N,1);

% simulate from age 6 to 12 investment
for t=6:T+1
    
    if (quo)
        % save pbar before price change
        pbar_age(:,t-1)=pbar;
        
        % save shares for t-1        
        if (~Cobb_Douglas)
            for m=1:2
                for e=1:2
                    share_model_age(m,e,t-1,:)=share_model(m,e,:);
                end
            end
        end
    end
    
    % logPsi_t
    logPsi=delta2*logPsi+qd_logPsi; % logPsi_{t}=delta2*logPsi_{t-1}+delta1*logX_{t-1}+logtheta
    
    % update X_t    
    if (t<=T)
        age(1:N)=t;
        nkid_0_5=(age<=5);
        
        % update regressors for the share parameters
        set_regressors;
        evaluate;
    end
end

if (quo)
    dlmwrite(['../temp/pbar_age' tech '.txt'],reshape(pbar_age,[N*T 1]));    
    
    if (~Cobb_Douglas)
        dlmwrite('../temp/share_model_age.txt', ...
                 reshape(share_model_age,[2*2*T*3 1]));
    end
end

for m=1:2
    
    j=max_j;
    mp=m-1;
    
    % log achievement at age 13
    j=j+1; % 10
    eval(['outcome(j,m)=mean(logPsi(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1; % 11
    aux=alpha.*logPsi; 
    eval(['outcome(j,m)=mean(aux(married==' int2str(mp) ...
          '),''omitnan'');']);    
        
end

if (quo)
    csvwrite(['../temp/quo' tech '.csv'],outcome)
else
    
    % compare with quo outcome
    outcome0=csvread(['../temp/quo' tech '.csv']);
    change=outcome-outcome0; % change
    change(1:7,:)=change(1:7,:)./outcome0(1:7,:)*100; % percentage change
    
    change(8,:)=change(8,:)*100; % 100*log change
    change(10,:)=change(10,:)*100; % 100*log change
    
    % consumption units: age 6 skill, age 5 investment
    change(9,:)=(exp( beta^(T-5+1)*delta2^(T-5)*change(9,:) )-1)*100;
    
    % consumption units: age 13 skill, ages 5-12 investment
    change(11,:)=(exp( (1-beta)/(beta^(5-1-T)-1)*change(11,:) )-1)*100;
    
    % save level and change
    csvwrite(['../temp/' name '_' num2str(pct_reduction) tech '.csv'],outcome)
    csvwrite(['../temp/change_' name '_' num2str(pct_reduction) tech '.csv'],change)
end
